Density matrix renormalization group approach of the spin-boson model 
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We propose a density matrix renormalization group approach to tackle a two-state system coupled 
to a bosonic bath with continuous spectrum. In this approach, the optimized phonon scheme is 
applied to several hundred phonon modes which are divided linearly among the spectra. Although 
DMRG cannot resolve very small energy scales, the delocalized-localized transition points of the 
two-state system are extracted by the extrapolation of the flow diagram results. The phase diagram 
is compared with the numerical renormalization group results and shows good agreement in both 
Ohmic and sub-Ohmic cases. 
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I. INTRODUCTION 

The density matrix renormalization group (DMRG) is 
an important tool for studying the strongly correlated 
systems in low dimensions ji*^ In the past decade, one sig- 
nificant limitation of DMRG — finite basis requirement in 
the model which involves infinite degree of freedom, e.g., 
phonon states, was circumvented by a controlled trunca- 
tion techniqueJi This technique is applied to many mod- 
els, typically, such as ID Holstein model,— ID Holstein- 
Hubbard modelf^ spin Peierls model, ^ and spin-boson 
modelf2. The main idea of this technique, controlled trun- 
cation, is realized by the density matrix approach which 
is useful for finding the most probable states of the trun- 
cated system. By its light, the infinite Hilbert space 
can be reduced to governable dimensions without sig- 
nificant loss of accuracy. However, the truncation tech- 
nique is originally designed for the systems involving just 
one phonon mode, i.e., the Einstein model. The direct 
application of this technique to the spin-boson model 
with a continuous spectrum of phonon modes, is not very 
successful.^ For instance, in the case of Ref. |6|, the num- 
ber of phonon modes were limited to = 18, the physics 
of this highly discrete model may be unreliable. Further- 
more, the number of states of each phonon mode kept 
is m = 2, the truncation error is relatively large and no 
convincing result on the delocalized-localized transition 
was found.' These limitation implies that, to handle the 
system with many phonon modes in a DMRG treatment, 
one needs to develop an improved truncation technique. 
This is the motivation of the present paper. 

Here, let us briefiy introduce the spin-boson model. 
The spin-boson model is an important toy model in the 
study of dissipative quantum systems. Its Hamiltonian 
is given by (set h= 



H 



A e 



'^ujia\ai + (Jz'^\i{ai + a\), (1) 



where the Pauli matrices ax,z describe a two-state sys- 
tem, a\ and are phonon creation and annihilation op- 
erators with frequencies oJi for the i-th phonon modes, e 
is an additional bias (asymmetry), A is the bare tunnel- 
ing splitting, and Ai represents the coupling between the 



two-state system and the i-th phonon mode. Generally, 
the so-called bath spectral function 



J{uj) = TT A^(5(cj — UJi) 



(2) 



completely determine the solutions of the spin-boson 
model. With an energy cutoff ojc, i.e., discards the high 
energy modes, the bath spectral function has a power-law 
form 



J(^) 



(3) 



where a is a dimensionless coupling constant which char- 
acterizes the dissipation strength, 0<s<l,s = l, and 
s > 1 represent sub-Ohmic, Ohmic, and super-Ohmic dis- 
sipation, respectively. The primary purpose of the spin- 
boson model is to study the effect of the environment 
on quantum tunneling of the two-state system. Here the 
environment is modelled as a collection of harmonic os- 
cillators, which serves as the origin of dissipation.^''* In- 
tuitively, the presence of the environment will make the 
tunneling particle as a "dressed" one, just like the elec- 
tron in a polaron-phonon (or exciton-phonon) system, 
and therefore its quantum tunneling decreases as the 
coupling increases. One important issue in spin-boson 
model is to study the phonon-induced localization (also 
stated as delocalized-localized transition), i.e., how quan- 
tum tunneling dies out as the coupling (or the dissipa- 
tion) increases^i^iSiiaiiiii^ii^iiiii^ii^ Such a delocalized- 
localized transition at T = is now considered as some 
kind of quantum phase transition called boundary phase 
transitioniiii^ 

In general, the Hamiltonian of the spin-boson model 
cannot be solved exactly, especially in the sub-Ohmic 
case.® The delocalized-localized transition has been 
widely studied by various methods with different approx- 
imations, yet a consensus on the delocalized-localized 
transition in sub-Ohmic case is still lacking. By inte- 
grating out the bath degrees of freedom, the spin-boson 
model was mapped to an Ising model and the localized 
transition was predicted to exist for s < 1 (i.e., in both 
Ohmic and sub-Ohmic cases). However, the path- 
integration with the so-called noninteracting blip ap- 
proximation (NIAP) and the adiabatic renormalization 
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predicted that no delocalized-localized transition hap- 
pens in the sub-Ohmic case*^ On the other hand, in the 
sub-Ohmic case, variational calculations, flow equation 
method, and other perturbation calculations predicted 
a discontinuous delocalized-localized transition^ii^ii^ii^ 
while the non-perturbative numerical renormalization 
group (NRG) calculation shows a continuous one.^^ Re- 
cently, the authors showed that the discontinuous tran- 
sition in the sub-Ohmic case obtained by the variational 
calculation is simply an artifact of the variational scheme 
due to the fact that the energy of the variational ground 
state can no longer be lower than the energy of the trial 
ground state (displaced-oscillator state) While this re- 
sult sheds some lights on the problem, the discrepancy 
between various treatments mentioned above has not 
yet been resolved. In addition, although the NRG ap- 
proach is regarded as the most powerful tool for treating 
the phase transition, the error due to discretization in 
numeric calculation has to be considered.^* Under this 
sense, the delocalized-localized transition is necessary to 
study by another non-perturbative method, i.e., DMRG. 
We hope that this paper will be helpful for resolving the 
discrepancy. 

The organization of this paper is as follows. In the 
following section, we propose a finite system DMRG ap- 
proach with controlled truncation technique to the spin- 
boson model, a thousand phonon modes can be treated. 
In Sec, mil we determine the DMRG parameters and dis- 
cuss the very small energy scales limitation of our treat- 
ment. Sec. IIVI suggests a extrapolation scheme to cir- 
cumvent the very small energy scales limitation. The 
delocalized-localized transition points of the spin-boson 
model, which is associated to the very small energy 
phonon modes, are obtained by extrapolation of "pseudo- 
critical" points. Conclusion is given in the last section. 



II. THE FINITE SYSTEM DMRG ALGORITHM 
OF THE SPIN-BOSON MODEL 

Here we present a finite system DMRG algorithm with 
the optimized truncation of multi-modes phonon space 
to treat the spin-boson model whose bosonic bath in- 
volves several hundred phonon modes. The key strategy 
of the algorithm is that one represents a single phonon 
mode as a site. The spin-boson model therefore be- 
comes a finite-size chain. In this case, finite system 
DMRG is naturally applied to this model since it is ap- 
propriate to reduce the environment error with sweep- 
ing processes.^" To reach this, we must divide the fre- 
quency spectrum into N intervals, i.e., [vi^i^i^i], where 
I = 1, . . . , TV, i/^ - Vi^i = Ui+i - Vi, 1^0=0, i>N ^ ujc^ 1, 
and LOi = {vi + Vi^i)/2. In other words, the frequency 
spectrum is divided linearly. The corresponding coupling 
parameters \i can be obtained by the spectral function 



© and ^ 

Xl = i r J{u)du = ^^^(-r^ - -f+i^). (4) 

TTJi/.^i 2(S-|-1) 

Suppose that Nb bare phonon states (|0), |1), . . . , |A^f,— 1)) 
are sufficient to represent one phonon mode accurately, 
therefore we can limit Nf, bare phonon states in each 
phonon site. Using the controlled reduction technique)^ 
the dimension of each phonon mode can be further re- 
duced to m where m < Nh. However, even though m, — 2 
is quite large for a dozen phonon modes, as in the treat- 
ment by Nishiyamai^ In this case, the number of phonon 
modes that can be treated is seriously restricted. Our so- 
lution to this problem is to truncate a set of phonon sites 
with the density matrix approach within each DMRG 
step, i.e., we do not optimized phonon modes individu- 
ally. The truncated multi-phonon sites can be contin- 
uously optimized by the sweeping of the finite system 
DMRG algorithm. Similar treatment was done by Fried- 
man in the study of spin-Peierls model. ^ With these pre- 
requisites, the finite system DMRG algorithm can be im- 
plemented in the following way. 

As the standard finite system DMRG algorithm which 
is used in Heisenberg model, ^ the first step of the algo- 
rithm is "warmup" . We must generate a series of phonon 
blocks for the subsequent sweeping processes of the finite 
system DMRG algorithm. For simplicity, we assume the 
number of phonon modes is odd and generate the blocks 
1 {N-l)/2 and {N+3)/2 N separately, where differ- 
ent numbers represent different phonon modes. With this 
simplification, all the phonon blocks 1 2, 1 3, . . . , 1 ^ 
{N-l)/2, {N+3)/2 ~ N, {N + 5)/2 - TV, . . . , A^- 1 ~ iV 
except for the blocks 1 and N which keep Nf, bare phonon 
states are limited to a 2 x 2 matrix because only the 
two-state system have been traced outj^ see Fig. [IJa). 
During the course of warmup, each phonon mode with 
Nf, bare phonon states is added to the preceding block. 
After a truncation with density matrix approach, one 
new block is generated. We shall show that Af, = 10 
is sufficient for the implementation of our algorithm in 
most cases. Note that, every block generated within the 
warmup processes must be stored in memory for later 
use. 

Secondly, the finite system DMRG algorithm is imple- 
mented as in Fig. [Ijb). The finite system DMRG algo- 
rithm is more or less the same as the standard algorithmii 
The main difference between the two algorithms is that 
we add one site within each DMRG step instead of two 
sites. It is because there are no interactions between the 
phonon blocks, the implementation of our algorithm is 
identical to the standard algorithm, and no further cor- 
rection is neededi^ For convenience, the two-state sys- 
tem can be placed on leftmost side or rightmost side 
on the chain, as to calculate the reduced density ma- 
trix and apply the traditional wave function transforma- 
tions technique to the finite system algorithm. ■^■^ Within 
each DMRG step, a phonon mode with Nb bare states is 
added. This new site is used to generate a new phonon 



3 



(a) 



^ \/ 
bare states \L 



optimized states 
J__2_ '^^1~2^_3_ 



(b) 



two-state system 



y 



1~(N-1)/ 2 



(N+1)/2 



C 



1 ~ (N-1)/2 J (8)^(8> r (N+3)/2~N ) 
(N+3)/2 



( 1~(N+1)/2 ) = (8> [ (N+5)/2~N ) (g) ^ 



1 ~ N-2 



n~)'8'=®r 



3 ~ N 



(N+1)/2 



( 1~(N-1)/2 )0 = (8>[ (N+3)/2~N ] ^ 



more than M = 20 30. In general, there are only 7 — 8 
largest eigenvalues in the reduced density matrix of the 
phonon block have significant values. The dependence 
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FIG. 1; (a) the warmup procedure of the finite system DRMG 
algorithm, where the numbers represent the phonon modes. 
This figure shows the warmup procedure of the left part 
phonon modes 1 ~ (TV — l)/2; see Fig. [ijb). The right 
part phonon modes can be obtained by similar fashion, (b) 
Systematic illustration of the finite system DMRG algorithm. 
This figure shows one sweep in the algorithm. 

block or optimize the old phonon block with M opti- 
mized states. Finally, the energies of the target states 
will converge after one or two sweeps are preformed. 
In summary, the algorithm can be proceeded as follows: 

1. warmup, generating a series of phonon blocks for 
subsequent sweeping processes; 

2. starting at the center phonon mode {N + l)/2, 
adding one phonon mode with A^f, bare states to 
the chain; 

3. performing the sweeping process to the whole 
chain; 

4. if the energies of the target states are not converged 
after a sweeping, then return to step (2). 

III. DISCUSSION AND THE LIMITATION OF 
THE ALGORITHM 

One important issue of the algorithm is that how to 
choose the parameters N^, Af , and the number of 
sweeps Ns- Unlike the spin 1/2 Heisenberg chain, there 
are no interactions between the phonon blocks, the num- 
ber of states kept per block is not quite large. Therefore, 
it may be possible to treat a thousand phonon modes 
while the number of states kept per block is never needed 



FIG. 2: Dependence of the ground state energy on the pa- 
rameters N,Nf,,M, and Ns for e = 0, s = 0.6, a = 0.1, and 
A = 0.1. (a) dependence on M for fixed N,Nb, and A^^; (b) 
dependence on Nh for fixed A'^, M , and Ns\ (c) dependence on 
the number of phonon modes N for fixed TV;,, Af, and A'^s; (d) 
dependence on Na, the number of sweeps versus the ground 
state energy, for fixed Nf,,N, and M. 

of the ground state energy on the parameters N,Ni,,M, 
and Ns for s = 0.6, a = 0.1, and A = 0.1 is shown in 
Fig. [2] (targeting the ground state only, but the following 
conclusions are also true for targeting both the ground 
state and the first excited state). It is worth noting that 
even though M = 10, Nj, = 6, and Ng — 1 can give rather 
the same results. However, the number of phonon modes 
N will highly affect the results. This is also the main 
limitation of our DMRG strategy. 

As indicated in Refs. d, [l^, and [13, the very small 
energy phonon modes are important for revealing the 
critical phenomena, e.g., the delocalized-localized tran- 
sition of the two-state system. However, the strategy 
of our algorithm needs linear discretization of the spec- 
trum which can not resolve very small energy scale. If 
one tries to apply a logarithmic discretization which is 
used in NRG to the DMRG algorithm, the energy levels 
of the Hamiltonian emerge a staircase-like aspect when 
one deals with the very small energy phonon modes; see 
Fig. [3l The DMRG scheme fails in this situation because 
the target states cannot be determined. This difficultly 
stems from the truncation strategy of the DMRG, say, 
it iteratively calculates the lowest eigenstates for finding 
the most probable states of the decimated system. How- 
ever, in practice, it is harsh for the iterative diagonaliza- 
tion routines being used by DRMG, such as Lanczos and 
Davidson, to converge when the staircase-like energy lev- 
els occur. In fact, the staircase-like energy levels also ap- 
pear in NRG calculations. Nevertheless, the truncation 
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FIG. 3: This figure is a schematic scaled energy spectrum of 
the DMRG or NRG calculation when one deals with the very 
small energy phonon modes. Here n is the number of energy 
levels. The vertical coordinate (energy) is scaled by A'^ and 
A'^ in DMRG and NRG, respectively. 



of NAE are qualitatively different within two regime 
a < ac and a > ac- Therefore, we assume that there 
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scheme of NRG, which retains the lowest-lying states di- 
rectly, is simply to implement in this situation. In other 
words, the performance of the standard diagonlization 
routine used by NRG will not be affected by the "shape" 
of the spectrum while DMRG needs iterative diagonliza- 
tion routine which converges arduously. 

DMRG cannot resolve very small energy scales, this 
limitation is serious. It implies that the critical cou- 
pling ttc cannot be determined due to the energy levels 
cannot reach to a fixed point without very small energy 
phonons^^ and the spin-spin correlation function cannot 
be calculated in very small energy scales. Furthermore, 
the effective tunneling splitting A,. = (ax) also is not ad- 
equate to identify the critical coupling Uc because it fails 
to characterize the tunneling in equilibrium in the sub- 
Ohmic case.^^ Our DMRG calculations have the same 
conclusion, namely, ^ when a > ac in the sub- 
Ohmic case and A^ — > when a ~> Oc (note that ac is a 
function of s and A) in the Ohmic case (not presented). 
Moreover, the entanglement entropy method proposed by 
Ref. [2^ recently, which is used to determine the critical 
couplings and performs very well in NRG, is not work- 
ing as expected in DMRG when the very small energy 
information is lacking. 



IV. EXTRACTING THE CRITICAL POINTS BY 
EXTRAPOLATION 

Now, we seek to show that the critical couplings ac 
can be determined by extrapolating the pseudo-critical 
couplings a'c which are extracted in a "DMRG flow" to 
thermodynamic limit. Similar to the energies in log- 
arithmic discretization of NRG which are falling off as 
A~^,^- the energies are falling off as iV~^ in linear dis- 
cretization. Therefore, one can target the ground state 
and the first excited state and scale the energy gap 
AE = i^excitod — ^-ground as NAE and plot the flow dia- 
gram NAE versus N. As one can see in Fig. HI the flows 



FIG. 4: Dependence of the scaled energy different between 
ground state and excited state NAE on the number of phonon 
modes N for a < a^, ce ~ a^, and q > Qc. Parameters are 
e = 0, s = 0.5, A''6 = 10, M = 20, and iV, = 2. 

exist a function a'c{N) which separates the two regimes, 
where TV is relatively small in comparison with thermody- 
namic limit. In practice, a'c{N) can be easily determined 
by a bisection process of two couplings a < a'^{N) and 
a > a'c{N) with the slope of a line segment consist of the 
scaled energies of two points [iV — 1, -I- 1]. We conceive 
that the pseudo-critical coupling a'^{N) will converge to 
the critical coupling ac when A/ — > 00 since the flxed 
points are reached. 

Accordingly, the extrapolation of the a'dN) versus 
1 /N curve determines the critical coupling ac at the limit 
of 1/iV 0. Figure shows the best-fit of the pseudo- 
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FIG. 5: Dependence of the pseudo-critical couplings 
ac( A'') (circle) and the best-fit values(line) on the inverse num- 
ber of phonon modes Parameters are e — 0,s = 
0.5, Nb = 10, M = 20, and iV^ = 2. 
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it turns out that ac ~ 0.09933. In fact, the best-fit 
curves a'^{l/N) are somewhat different for sub-Ohmic 
and Ohmic dissipation. It is related to the fact that the 
transition in the sub-Ohmic case is characterized by a 
quantum critical fixed point in contrast to the Ohmic 
casej^ Since we cannot find a formula to fit all the cases, 
the extrapolations are done by the simplest polynomial 
fitting. 

Admittedly, it might be doubted that if N is small, 
a'^ could be inappropriate for extrapolation since it is 
inconsistent with the a'^ which are obtained with large 
N. For instance, there are some cases show that the 
infinite system DMRG is a better choice to tackle this 
problem However, on the one hand, our strategy 
of the DMRG in the spin-boson model limits the im- 
plementation of the algorithm. In order to "insert" the 
spin-boson model into the DMRG algorithm, one must 
cut the spectrum of the bosonic bath to finite number 
of pieces and therefore the spin-boson model becomes a 
finite-size chain. Before performing the linear discretiza- 
tion, the number of sites N and the coupling constant 
Ai must be determined. Naturally, it brings about the 
finite system DMRG algorithm to handle this model and 
a "real" infinite DMRG algorithm is difficult to imple- 
ment in practice. On the other hand, in our finite system 
DMRG solution, the number of phonon modes treated is 
relatively large. We carefully check the calculations and 
find that when we calculate the a'^ with < 0.005, 
the a'^ are always monotonic. Hence, the extrapolations 
are safe and correct in our treatment. Furthermore, ref- 
erences [2^ and [2^ also performed an extrapolation of the 
number of states kept, but the result shown in Fig. [2] and 
the fact of non-interacting phonon blocks guarantee that 
this quantity is not significant in our calculation notwith- 
standing. 

Using the extrapolation scheme, the phase boundary 
for the delocalized-localized transition of the spin-boson 
model for e = and A = 0.1 is shown in FiglHl where the 
result by NRG is also shown for comparison. It can be 
found that the DMRG data are consistent with the NRG 
data quite well. It also shows that the NRG data are 
always larger than the DMRG data. Indeed, however, 
the results of NRG can be extrapolated to thermody- 
namics limit, namely, one takes the NRG discretization 
parameter A ^ 1, and smaller critical couplings can be 
obtainedfi^ In other words, both NRG and DMRG show 
that the critical couplings before extrapolating to ther- 
modynamics limit are always larger than the true criti- 
cal couplings. This implies that the error of discretiza- 



tion on determining the critical coupling is to lower the 
true ac, but not to heighten the ac as claimed in Ref. 
H In addition, the inset of Fig. [HI also assures that our 
DMRG calculations for Ohmic case are consistent with 
the NRG result^^ and the well-known renormalization 
group resultji i.e., ac = 1 + 0{A/lJc)- 
V. CONCLUSION 

At the end of this paper, we want to compare 
the ground state energy obtained by DMRG with the 
variational ground state energy obtained by displaced- 
oscillator approach. As we mentioned in the Introduc- 
tion, the variational ground state will fail at a cer- 
tain point since the variational ground state is not yet 
stable J^ii^ The comparison of the ground state energy is 
shown in Fig. [T] It shows that the DMRG ground state 
energies are always lower than that of the variational cal- 
culations. It is also clear that the DMRG ground state 
energy approaches to variational ground state energy for 
small a and displaced-oscillator ground state energy for 
large a, and no discontinuous effective tunneling splitting 
and ground state energy are observed by the DMRG cal- 
culations. 

In conclusion, we have proposed a finite system DMRG 
algorithm to deal with the spin-boson model. This algo- 
rithm is much more powerful than the preceding study 
of this topic because it can treat more than a thousand 
phonon modes.''^ In fact, we have tried to calculate lO'' 
and 10^ phonon modes. Unfortunately, our 32-bit system 
is unable to tackle phonon modes on the 10^ magnitude 
due to out of memory. This difficultly, of course, can be 
resolved in 64-bit system or storing the data in hard disk. 
Moreover, we obtain the critical couplings ac by extrap- 
olating the pseudo-critical couplings to thermodynamic 
limit. The phase diagram is compared with NRG and 
shows good agreement in both sub-Ohmic and Ohmic 
cases. 
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